home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / matrix_sub.cc < prev    next >
Encoding:
Text File  |  1995-12-18  |  6.2 KB  |  245 lines  |  [TEXT/ttxt]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Linear Algebra Package
  6.  *
  7.  *        Basic linear algebra operations, levels 1 & 2
  8.  *        Operations on a single row, column, or the diagonal
  9.  *                   of a matrix
  10.  *
  11.  * $Id: matrix_sub.cc,v 3.1 1995/01/31 17:27:47 oleg Exp oleg $
  12.  *
  13.  ************************************************************************
  14.  */
  15.  
  16. #include "LinAlg.h"
  17. #include <math.h>
  18.  
  19.  
  20. /*
  21.  *------------------------------------------------------------------------
  22.  *            Messing with a single column of a matrix
  23.  */
  24.  
  25.                  // For every element, do `elem OP value`
  26. #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                \
  27.                                     \
  28. void MatrixColumn::operator OP (const VALTYPE val)            \
  29. {                                    \
  30.   matrix.is_valid();                            \
  31.                                     \
  32.   register REAL * cp = ptr;        /* Column ptr */        \
  33.   while( cp < ptr + matrix.nrows )                    \
  34.     *cp++ OP val;                            \
  35. }                                    \
  36.  
  37. COMPUTED_VAL_ASSIGNMENT(=,REAL)
  38. COMPUTED_VAL_ASSIGNMENT(+=,double)
  39. COMPUTED_VAL_ASSIGNMENT(*=,double)
  40.  
  41. #undef COMPUTED_VAL_ASSIGNMENT
  42.  
  43.  
  44.                 // Assign a vector to a matrix col
  45. void MatrixColumn::operator = (const Vector & vec)
  46. {
  47.   matrix.is_valid();
  48.   vec.is_valid();
  49.  
  50.   if( matrix.row_lwb != vec.row_lwb || matrix.nrows != vec.nrows )
  51.     matrix.info(), vec.info(),
  52.    _error("Vector cannot be assigned to a column of the mentioned matrix");
  53.  
  54.   register REAL * cp = ptr;        // Column ptr
  55.   register REAL * vp = vec.elements;    // Vector ptr
  56.   while( cp < ptr + matrix.nrows )
  57.     *cp++ = *vp++;
  58.  
  59.   assert( vp == vec.elements + vec.nrows );
  60. }
  61.  
  62.                 // Assign a matrix column to a vector
  63. Vector& Vector::operator = (const MatrixColumn & mc)
  64. {
  65.   is_valid();
  66.   mc.matrix.is_valid();
  67.  
  68.   if( mc.matrix.q_row_lwb() != row_lwb || mc.matrix.q_nrows() != nrows )
  69.     mc.matrix.info(), info(),
  70.    _error("Can't assign a column of the matrix above to the vector");
  71.  
  72.   register REAL * cp = mc.ptr;        // Column ptr
  73.   register REAL * vp = elements;    // Vector ptr
  74.   while( vp < elements + nrows )
  75.     *vp++ = *cp++;
  76.  
  77.   assert( cp == mc.ptr + mc.matrix.q_nrows() );
  78.  
  79.   return *this;
  80. }
  81.  
  82.  
  83. /*
  84.  *------------------------------------------------------------------------
  85.  *             Messing with a single row of a matrix
  86.  *  Keep in mind the matrix elements are arranged in columns!
  87.  */
  88.  
  89.  
  90.                  // For every element, do `elem OP value`
  91. #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                \
  92.                                     \
  93. void MatrixRow::operator OP (const VALTYPE val)                \
  94. {                                    \
  95.   matrix.is_valid();                            \
  96.                                     \
  97.   register REAL * rp = ptr;        /* Row ptr */            \
  98.   for(; rp < ptr + matrix.nelems; rp += inc)                \
  99.     *rp OP val;                                \
  100. }                                    \
  101.  
  102. COMPUTED_VAL_ASSIGNMENT(=,REAL)
  103. COMPUTED_VAL_ASSIGNMENT(+=,double)
  104. COMPUTED_VAL_ASSIGNMENT(*=,double)
  105.  
  106. #undef COMPUTED_VAL_ASSIGNMENT
  107.  
  108.  
  109.                 // Assign a vector to a matrix row
  110.                 // The vector is considered row-vector
  111.                 // to allow the assignment in the strict
  112.                 // sense
  113. void MatrixRow::operator = (const Vector & vec)
  114. {
  115.   matrix.is_valid();
  116.   vec.is_valid();
  117.  
  118.   if( matrix.col_lwb != vec.row_lwb || matrix.ncols != vec.nrows )
  119.     matrix.info(), vec.info(),
  120.    _error("Transposed vector cannot be assigned to a row of the matrix");
  121.  
  122.   register REAL * rp = ptr;        // Row ptr
  123.   register REAL * vp = vec.elements;    // Vector ptr
  124.   for(; rp < ptr + matrix.nelems; rp += inc)
  125.     *rp = *vp++;
  126.  
  127.   assert( vp == vec.elements + vec.nrows );
  128. }
  129.  
  130.                 // Assign a matrix row to a vector
  131.                 // The matrix row is implicitly transposed
  132.                 // to allow the assignment in the strict
  133.                 // sense
  134. Vector& Vector::operator = (const MatrixRow & mr)
  135. {
  136.   is_valid();
  137.   mr.matrix.is_valid();
  138.  
  139.   if( mr.matrix.col_lwb != row_lwb || mr.matrix.ncols != nrows )
  140.     mr.matrix.info(), info(),
  141.    _error("Can't assign the transposed row of the matrix to the vector");
  142.  
  143.   register REAL * rp = mr.ptr;        // Row ptr
  144.   register REAL * vp = elements;    // Vector ptr
  145.   for(; vp < elements + nrows; rp += mr.inc)
  146.     *vp++ = *rp;
  147.  
  148.   assert( rp == mr.ptr + mr.matrix.nelems );
  149.  
  150.   return *this;
  151. }
  152.  
  153. /*
  154.  *------------------------------------------------------------------------
  155.  *            Messing with the matrix diagonal
  156.  */
  157.  
  158.                  // For every element, do `elem OP value`
  159. #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                \
  160.                                     \
  161. void MatrixDiag::operator OP (const VALTYPE val)            \
  162. {                                    \
  163.   matrix.is_valid();                            \
  164.                                     \
  165.   register REAL * dp = ptr;        /* Diag ptr */            \
  166.   register int i;                            \
  167.   for(i=0; i < ndiag; i++, dp += inc)                    \
  168.     *dp OP val;                                \
  169. }                                    \
  170.  
  171. COMPUTED_VAL_ASSIGNMENT(=,REAL)
  172. COMPUTED_VAL_ASSIGNMENT(+=,double)
  173. COMPUTED_VAL_ASSIGNMENT(*=,double)
  174.  
  175. #undef COMPUTED_VAL_ASSIGNMENT
  176.  
  177.  
  178.                 // Assign a vector to the matrix diagonal
  179. void MatrixDiag::operator = (const Vector & vec)
  180. {
  181.   matrix.is_valid();
  182.   vec.is_valid();
  183.  
  184.   if( ndiag != vec.nrows )
  185.     matrix.info(), vec.info(),
  186.    _error("Vector cannot be assigned to the matrix diagonal");
  187.  
  188.   register REAL * dp = ptr;        // Diag ptr
  189.   register REAL * vp = vec.elements;    // Vector ptr
  190.   for(; vp < vec.elements + vec.nrows; dp += inc)
  191.     *dp = *vp++;
  192.  
  193.   assert( dp < ptr + matrix.nelems + inc );
  194. }
  195.  
  196.  
  197.                 // Assign the matrix diagonal to a vector
  198. Vector& Vector::operator = (const MatrixDiag & md)
  199. {
  200.   is_valid();
  201.   md.matrix.is_valid();
  202.  
  203.   if( md.ndiag != nrows )
  204.     md.matrix.info(), info(),
  205.    _error("Can't assign the diagonal of the matrix to the vector");
  206.  
  207.   register REAL * dp = md.ptr;        // Diag ptr
  208.   register REAL * vp = elements;    // Vector ptr
  209.   for(; vp < elements + nrows; dp += md.inc)
  210.     *vp++ = *dp;
  211.  
  212.   assert( dp < md.ptr + md.matrix.nelems + md.inc );
  213.  
  214.   return *this;
  215. }
  216.  
  217. /*
  218.  *------------------------------------------------------------------------
  219.  *           Multiplications with the diagonal matrix
  220.  */
  221.  
  222.                 // Multiply a matrix by the diagonal
  223.                 // of another matrix
  224.                 // matrix(i,j) *= diag(j)
  225. Matrix& Matrix::operator *= (const MatrixDiag& diag)
  226. {
  227.   is_valid();
  228.   diag.matrix.is_valid();
  229.  
  230.   if( ncols != diag.ndiag )
  231.     info(), diag.matrix.info(),
  232.    _error("Matrix cannot be multiplied by the diagonal of the other matrix");
  233.  
  234.   register REAL * dp = diag.ptr;    // Diag ptr
  235.   register REAL * mp = elements;    // Matrix ptr
  236.   register int i;
  237.   for(; mp < elements + nelems; dp += diag.inc)
  238.     for(i=0; i<nrows; i++)
  239.       *mp++ *= *dp;
  240.  
  241.   assert( dp < diag.ptr + diag.matrix.nelems + diag.inc );
  242.  
  243.   return *this;
  244. }
  245.